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Abstract. We study by molecular dynamics the structural properties of single layer 
h-BN in comparison to graphene. We show that the Tersoff bond order potential 
developed for BN by Albe, Moller and Heinig gives a thermally stable hexagonal single 
layer with a bending constant k = 0.54 eV at T = 0. We find that the non-monotonic 
behaviour of the lattice parameter, the expansion of the interatomic distance and the 
growth of the bending rigidity with temperature are qualitatively similar to those of 
graphene. Conversely, the energetics of point defects is extremely different: instead 
of Stone Wales defects, the two lowest energy defects in h-BN involve either a broken 
bond or an out of plane displacement of a N atom to form a tetrahedron with three B 
atoms in the plane. We provide the formation energies and an estimate of the energy 
barriers. 

PACS numbers: 65.40.De, 65.80.-g, 68.35.Dv, 61.72.Bb 
1. Introduction 

The interest in two-dimensional (2D) crystals is rapidly extending to other materials 
than graphene, often combined to form man-made heterostructurcs[lJ. Hexagonal 
boron-nitride (/i-BN), which has similar structural properties to graphene but is an 
insulator with a gap of ~ 5 — 6 eVp], is one of the promising dielectric materials for 
integration in hybrid graphene devices Deposition of graphene on /i-BN is found to 
improve the transport properties of graphene possibly due to suppression of scattering 
from out of plane ripples P]. The growth of thin BN films is known to require a 
certain amount of ion irradiation which in turn leads to the detrimental formation 
of point defects |5]. The presence of defects lead to special features in x-ray core level 
spectroscopy [51 E] which have stimulated the study of selected defects in /i-BN by first- 
principles [3 El |9l [10]. Another interest in the defect structure of /i-BN is the possibility 
to create suitable chemisorption sites for molecules as proposed for graphene [TT]. In 
this paper, we address the structural stability, thermal expansion and defect formation 
in /i-BN by means of Molecular Dynamics (MD) based on a classical description of 
interatomic interactions. The success of this approach for graphene, particularly for 
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temperature dependent properties, has been due to the existence of several accurate 
models of interactions [T2| [T3| HU [15]. These so called bond-order potentials for carbon, 
pioneered by Tersoff[T2], have the important feature of being reactive, namely to allow 
change of coordination. For BN, Albe, Moller and Heinig[T6] have developed a Tersoff 
potential that is supposed to describe well the bulk BN phases and the defect formation 
energy and has been used to study the effect of irradiation in single layer /i-BN[T7]. To 
our knowledge, however, no comprehensive study of the structural stability and defect 
formation energies based on this approach exists up to date. The goal of this article 
is therefore two-fold: on the one hand, to study the structure of single layer /i-BN in 
comparison to graphene and report an extensive search for low energy point defects 
and, on the other hand, provide a wide set of results that can be later tested against 
experiments or more sophisticated calculations. Validation of this potential is necessary 
also in view of a possible extension to a reliable potential for hybrid systems formed by 
B, N and carbon. We find that the thermal expansion and bending rigidity of /i-BN 
behave similarly to graphene whereas the type and formation energy of lattice defects 
is very different. The paper is organized as follows: in section [2] we describe our method 
and studied samples, in section [3] and |4] we present our results for the temperature 
dependence of thermal expansion and bending rigidity of crystalline single layer /i-BN 
respectively. In section [5] we show how, by raising the temperature up to melting, 
we identify several possible lattice defects for which we give the formation energy. In 
section |6] we give a summary and perspectives of our work. 

2. Method 

We have performed MD simulations using the lammps [TH] software package and the 
BN interatomic potential of Albe, Moller and Heinig [16]. This potential gives at T = 
the /i-BN equihbrium lattice parameter a = a/S-Rq =2.532 A |16] where Rq is the B- 
N interatomic distance. Experimentally, the lattice parameter of /i-BN is 2.504 A and 
decreases with increasing temperature [19]. This means that the potential overestimates 
this value, resulting in a larger mismatch of 2.9% with the lattice parameter of 2.46 A 
for graphene, instead of the 1.8% experimental value. 

To study the thermal stability and lattice expansion we have performed simulations 
starting with a flat sheet of /i-BN consisting of n = 9800 atoms (sample A) corresponding 
to sample sides = 177.3 A and Ly = 153.5 A and periodic boundary conditions. 
These simulations were performed in the isobaric-isothermal (nPT) ensemble using a 
Berendsen thermostat with damping parameter = 0.1 ps to control the temperature 
and a Berendsen barostat [20j to impose an external pressure P = 0. A smaller sample of 
n = 1296 atoms {L^ = 45.59 A and Ly = 78.96 A) (sample B) was used to study defects 
by performing simulations in the canonical (nVT) ensemble using the T = lattice 
parameter. We have found that simulations with the Nose-Hoover thermostat yielded 
similar values for the lattice expansion but had undesired fluctuations in the enthalpy. 
The time step of all simulations was set to At = 0.1 fs and integration was performed 
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using the standard velo city- Ver let algorithm. For each simulation with sample A, a total 
of 2.5 millions time steps were taken, where the last 2.0 millions (200 ps) were used for 
averaging. 

3. Lattice Parameter 

One of the interesting quantities is the lattice thermal expansion. For graphene, 
simulations based on the bond order potential LCBOPII [21] have found that the lattice 
parameter first decreases with increasing temperature up to ~900 K and then increases 
at higher temperature. Calculations based on the quasi-harmonic approximation predict 
instead a contraction of the lattice parameter at least up to 2500 Kj23j. Experimentally, 
data up to 400 K for graphene [23] and up to room temperature for /i-BN [19] confirm 
a contraction. We distinguish three different structural parameters, namely the in- 
plane lattice constant a determined from the equilibrium box size obtained at constant 
pressure P = 0, the B-N nearest neighbour distance R and the nearest neighbour 
distance projected in the xy-plane R^y At T = K, for a completely fiat sheet of 
/i-BN, R = Rxy = a/\^. In figure l]we show the calculated temperature dependence 
of these three quantities for single layer /i-BN. We find that the BN nearest neighbor 
distance increases linearly up to about 1500 K with a slope of 1.2 ■ 10~^ A K~^, which is 
almost twice as large as the value 6.5-10"^ A found in [22] for freestanding graphene 
using ab initio MD simulations. For the lattice parameter we find thermal contraction 
up to 1200 K, similarly to prediction for graphene[2I]. For higher temperatures, the 
lattice parameter grows less rapidly than for graphene and, in the studied temperature 
range up to 2500 K, it never gets larger than the value at T = K. At T 7^ 0, the lattice 
parameter a is not equal to y/SR^y due to fluctuations in the out of plane direction. 

4. Bending rigidity k, 

We analyse the height fluctuations h of the single layer /i-BN crystal by using the theory 
of membranes in the continuum (see for instance [23 |26l [23 [2H1 [22] )• Key quantities in 
this theory are the correlation functions of the height fluctuations and of the normals. 
The normal-normal correlation function G{q) in the harmonic approximation is given 



where S = L^Ly/n is the area per atom, k is the bending rigidity and the subscript zero 
indicates that averages are taken in the harmonic approximation, namely neglecting 
the coupling of in-plane and out of plane fluctuations in the stress tensor [26]. In the 
limit of slowly varying height fluctuations, one can show that G{q) is related to the 
height-height correlation function H{q) by: 




(1) 




(2) 
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yielding, in the harmonic approximation, 

(3, 

It has been shown [28 j that Equation ^ is well reproduced when using numerical 
results from atomistic simulations only if the height fluctuations are calculated by 
averaging over the nearest neighbours heights: 

h=l(ho + hh^ + hf} + h^)] (4) 



2 V 3 

where ho is the z coordinate of one atom and hi the z coordinates of its three nearest 
neighbours. In the top panel of figure [2| we compare G{q) to q^H[q) calculated by MD 
at T = 300 K. We notice that these functions indeed coincide for g < 1 A^^. Above 
this value the continuum approximation used in the theory of membranes breaks down 
and deviations from power-law behaviour occur, resulting in a peak at the position of 
the Bragg peak q = 47r/\/3a = 2.86 A"^. 

The theory of membranes [26| 129] predicts deviations from harmonic behaviour for 
wavevectors smaller than 



where Y is the bulk modulus. In the longwavelength limit, for q < q*, the correlation 
functions bend to a lower exponent. This feature reduces the divergence of out of plane 
fluctuations and stabilizes the 2D crystal [221 ED]- In the top panel of figure |2] we see 
that deviations of the calculated points from Go{q) start at q* ~ 0.24 A~^ as found for 
graphene at the same temperature [2HI EI], suggesting a similar ratio Y/k. 

In the bottom panel of figure [2] we show H(q)/n calculated by MD at different 
temperatures. The bending rigidity k can be obtained by a fit of the slope of 
the calculated H{q) to Equation ([s]) in the range of q vectors where the harmonic 
approximation applies, yielding the temperature behaviour shown in figure |3} The 
value of K at T = can be found by direct total energy calculations of nanotubes and 
extrapolation to the limit of infinite radii[32j. This procedure yields n = 0.54 eV, a 
value lower than the value k = 0.82 eV for graphene[27], meaning that /i-BN is easier 
to bend. 



5. Defects 



Classical MD simulations are very suitable to search for possible distortions of the 
lattice. By raising the temperature we have observed the formation of some defects that 
arise in the melting process which occurs spontaneously at T ~ 4000 K. We have then 
studied the energetics of these defects and of others defects that are known to occur in 
graphene and in semiconductors. We distinguish two kinds of point defects, those where 
the number of atoms remains the same and those where one or more atoms are removed 
(vacancies). For the first group of defects we can calculate the formation energy as: 

Ep = E defect — Eperfect, (6) 
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where Edefect and Eperfect are the total energies of the sample with and without the 
defect. The cohesive energy is E^oh = Eperfect/n = —6.4166 eV. For these calculations 
we use a sample of n = 1296 atoms (sample B) which is cooled down to T = K in 
the nVT ensemble. In these calculations all defects were artificially created, thus not 
always observed in high temperature simulations. 

First of all we consider Stone-Wales (SW) defects |33] where a 90 degrees rotation 
of a pair of atoms transforms four hexagons into two pentagons and two heptagons (see 
figure |4^)). In graphene SW defects are known to have the lowest formation energy ~ 
4.7 eV [31] • Although SW defects have been theoretically studied in both nanotubes and 
single layer /i-BN[7l[H], they have not been found experimentally [55] . Also in our melting 
simulations we have not observed their formation. Recently, it has been suggested [7] 
that SW defects in graphene are further stabilized by a sine-like or cosine-like buckle 
and that this finding should hold also for other hexagonal single-layer crystals. By 
starting with a flat layer with a SW defect, and raising the temperature to 10 K we find 
a spontaneous buckling to one of the two structures SWi and SW2 shown in figure ^p) 
and figure [Ij^) respectively which resemble the sine-like buckle proposed in [7]. We found 
the cosine-like buckle to be unstable and deform to the structure SW2- The formation 
energies (table [T]) of SWi is slightly lower and both are almost twice the value of a SW 
in graphene [H [51] . 

While raising the temperature close to melting, we observe first two defects that 
have no counterpart in graphene, and turn out to have the lowest formation energy. In 
figure |5^) we show the defect that we call BB-defect because it results from a broken 
BN bond. The B atom moves to a metastable state with a slightly longer BN bond 
length (1.53 A instead of 1.47 A) with the two neighbours and the formation of two 
loose B-B bonds with a large bond length (1.9 A). No extended deviations in the z- 
direction occur for this defect (see table [T]). The formation energy (4.4 eV) is only 
0.1 eV higher than the tetrahedron defect (figure ^)) where a N atom pops up and 
form a NB3 tetrahedron with the three B atoms in the plane. Also this defect does 
not lead to out of plane deformations in the surrounding. It has been suggested that a 
BN3 tetrahedron structure with one N atom on top could explain some features in core 
level spectroscopy [5] E]- We find that this structure is unstable and that the B atom 
abandons the layer if lifted from the plane. In figure |6] we show the energy profile for 
the BB and tetrahedron defects calculated by the nudged elastic band (NEB) method, 
as implemented in lammps[T8|. We see that the barrier for the BB defect is 2.9 eV 
higher than that of the tetrahedron defect. In figure [5]d) we show also the antisite defect 
which is commonly found in semiconductors. As shown in table [T] its formation energy 
is intermediate between those of the BB and tetrahedron defects and those of the SW 
defects. 

In table [T] we also give the maximum out of plane displacement and the average 
value of out of plane displacements in the whole sample. We can see that, with the 
exception of the antisite defect, the formation energy seems to be related to the out of 
plane distortions, as suggested for grain boundaries in graphene [36]. 
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Lastly, we consider the various vacancies occurring when one or more atoms are 
removed. We define the vacancy giving as subscript the removed atoms, e.g. V^r as a 
system missing one N atom. Apart from the Vbn vacancy that keeps the stoichiometry, 
for these cases the evaluation of the formation energy requires to know the chemical 
potential of bulk phases of the constituents [37]. Different formation energies can result 
from use of different reference systems. Moreover, the chemical potential is often 
approximated by the cohesive energy, i.e. its value at T = 0|37]. For BN, the formation 
energy of the Vbn, V^b+n and Vb+sn vacancies has been calculated ab-initio in p3] by 
use of the chemical potentials of bulk metallic boron and solid nitrogen. As discussed 
in Ref . |38] . this approach is not reliable for a phenomenological potential like the one 
we are using. Therefore in table |2] we just give the energy difference AE between the 
energy of the perfect sample {Eperject = nEcoh) and the one of the sample relaxed after 
creation of the vacancy. We notice that Vb and Vn have the same AE because they 
both involve three BN broken bonds. Once the vacancy is created, no new bonds are 
formed. The three atoms surrounding the vacancy remain two-fold coordinated with 
the same bond angle. The bond length with the two nearest neighbours contracts from 
1.46 A to 1.44 A. 

Since, the single vacancies Vn and Vb have the same AE, we have a qualitative 
indication of the formation energy by subtracting from AE the cohesive energy for each 
removed atom, irrespective of its nature, namely A = AE — n^Ecoh, where is the 
number of missing atoms. 

Also for Vbat, Vsb+n and Vb+sn we find that no new bonds are formed after creation 
of the vacancy and that the structural changes are neglig ible. In table [2] we give the 
corresponding values of AE and A. To establish whether these are indeed the lowest 
energy structures, we have brought the atoms around the vacancy closer to each other 
inducing the formation of new bonds re-establishing three-fold coordination as shown in 
figure [7]for Vbn, V^b+n and Vb+sn- All these structures remain bonded after relaxation 
but have higher energy than the ones with no rebonding. This is due to the strong N-N 
bonds of only ^ 1.01 A in Vbat and V^b+n that cause strong local out of plane distortions 
and hinder the ring formation. Indeed, the ring structure of the V^b+n vacancy shown 
in figure [7]d) is obtained by construction and does not occur spontaneously whereas the 
Vb+3n, with no N-N bonds, may form spontaneously at finite temperature since it has 
an energy only marginally higher than the one without rebonding. 

6. Summary and conclusions 

In summary, we have presented an extended study of the structural properties of single 
layer /i-BN by means of MD simulations based on the interatomic potential developed by 
Albe, MoUer and Heinig and compared our findings to those for graphene. Validation 
of the results given by this potential opens the way to the development of a reliable 
potential capable to deal with hybrid BN-graphene structures. We find that the non- 
monotonic behaviour of the lattice parameter, the expansion of the interatomic distance 
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and the growth of the bending rigidity with temperature are quahtatively similar to those 
of graphene. Conversely, the energetics of point defects is extremely different: Stone 
Wales defects have formation energy twice as large as the two lowest energy defects in 
/i-BN which involve either a broken bond or an out of plane displacement of a N atom to 
form a tetrahedron with three B atoms in the plane. We have also studied the antisite 
defect and vacancies formed by one, two and four atoms. We hope that our results will 
stimulate further research on this topic. 
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Table 1. Formation energy Ep = Edefect — Eperfect of defects that do not change 
the number of atoms. The height difference h along the z-axis between the highest and 
lowest points is also given. In brackets we give the standard deviation [{z^) — (z) ]^/^ 



evaluated 


over the whole sample. 


Defect 


Ep (eV) 


h{k) 


SWi 


8.6 


2.26 (0.22) 




8.8 


2.04 (0.28) 


BB 


4.4 


1.12 (0.13) 


Tetra 


4.3 


1.43 (0.20) 


Anti-site 


6.3 


1.41 (0.10) 



Table 2. The number riy of atoms removed to create the vacancy and the energy 
difference Ai? between the perfect sample and the sample with a defect for various 
vacancies. The quantity A = AE — riyEcoh is given for the minimal energy structure 
without rebonding (see text), while A' is the same quantity for the defects after ring 
formation (see figure [?]) . For the latter defects we also give the height difference h 
along the z-axis between the highest and lowest points and in brackets the standard 
deviation [(z^) — {z) ]^^^ evaluated over the whole sample. 



Defect 




AE (eV) 


A (eV) 


A' (eV) 


h{A) 


Vb 


1 


11.7 


5.3 






Vn 


1 


11.7 


5.3 






Vb+n 


2 


19.7 


6.9 


8.8 


2.06 (0.35) 


Vb+3N 


4 


36.0 


10.4 


10.9 


1.42 (0.31) 


VsB+N 


4 


36.2 


10.6 


16.1 


2.63 (0.22) 
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1.50 




Figure 1. Temperature dependence of lattice parameters calculated by MD in the 
nPT ensemble for sample A consisting of n = 9800 atoms. (Top panel) Average BN 
interatomic distance R. A linear fit to the calculated points with R ^ Rq + aT where 
Rq =1.462 A and a = 1.2 • lO"'^ A K'^ describes well the results up to ~ 1500K. 
(Bottom panel) Lattice parameter a (left y-axis) and BN interatomic distance projected 
in the xy-plane Rxy (right y-axis). Notice that left and right y-axis differ by a factor 

Vs. 
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Figure 2. (Top panel) The two correlation functions G{q)/n Equation (III and 
q^H{q)/n Equation ([3| at T = 300K coincide as expected from Equation ^2]) for 
wave vectors (? < 1 yielding the same value of n from a best fit to the harmonic 
approximation Equation ^ (dotted line) in the range q — 0.5 — 0.9 A^^. Deviations 
from power-law behaviour occur for q }^ 1 A^^ close to the Bragg peak position 
at q = An/V^a = 2.86 A~^ In the long wavelength hmit, q < q* « 0.24 A~\ the 
correlation functions have a power-law behaviour with a smaller exponent |31). (Bottom 
panel) H{q)/n for T = 300 K, T = 1500 K and T = 2500 K. 



> 1.0 

y 0.9 




500 1000 1500 2000 2500 

Temperature (K) 



Figure 3. Bending rigidity function of temperature. The T = K value 

K = 0.54 eV is found by direct total energy calculations (see text). 
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c) 



Figure 4. Stone- Wales defects in BN. A pair of atoms is rotated by 90 degrees to 
form two pentagons and two heptagons (a). Side view of the two buckled structures 
with similar formation energy (see text): b) SWi . The bond lengths are 1) 1.47 A, 
2) 1.71 A, 3) 1.75 A and 4) 1.58 A. c) SW2. The bond lengths 1 and 3 change slightly 
with respect to SWI: 1) 1.45 A and 3 1.72 A. 




Figure 5. top and side view of a) the BB defect, b) the antisite defect and c) the 
tetrahedron defect. Selected values of the interatomic distances: 1) 1.53 A 2) 1.90 A 
3) 1.60 A 4) 1.47 A 5) 1.59 A 6) 1.61 A 7) 1.41 A. In all three cases the defects only 
cause a local distortion in the z-direction leaving the rest of the sample relatively flat 
(see table [1]). 
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Figure 6. The energy profile calculated using the nudged elastic band method for 
the BB and the tetrahedron defect. The reaction coordinate is only well defined at 
the points where it represents the perfect sample and 1 for the final sample with the 
defect. The formation energy is defined to be Ep — Edefect — Eperfect- The energy 
difference between the two highest points is 2.9 eV. 
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Figure 7. Top and side view of a) Vbn b) Vsb+at c) Vb+zn where rebonding and 
ring formation is obtained by construction (see text and table[2]). In both a) and b) 
short N-N bonds are created leading to strong out of plane distortions. Selected values 
of interatomic distances: 1) 1.01 A 2) 1.91 A 3) 1.92 A 4) 1.47 A 5) 1.90 A 6) 1.49 A 
7) 1.53 A 8) 1.58 A 9) 1.55 A 10) 1.52 A. 



